Consignas

A.- Para la base de datos seleccionada genere una muestra aleatoria estratificada y balanceada por variedad de vino de tamaño n = 2000 utilizando como semilla los últimos tres dígitos del DNI/PASAPORTE.

B. De ahora en más trabaje con esta base de datos para el resto del trabajo práctico. Realice los procedimientos que se detallan a continuación acompañando los procedimientos de los gráficos que considere adecuados.

datos_tp <-read_excel("DatosTP1.xlsx") %>% mutate_at('variedad', as.factor)
glimpse(datos_tp)
## Rows: 6,497
## Columns: 13
## $ `acidez fija`               <dbl> 7.0, 6.3, 8.1, 7.2, 7.2, 8.1, 6.2, 7.0, 6.…
## $ `acidez volátil`            <dbl> 0.27, 0.30, 0.28, 0.23, 0.23, 0.28, 0.32, …
## $ `ácido cítrico`             <dbl> 0.36, 0.34, 0.40, 0.32, 0.32, 0.40, 0.16, …
## $ `azúcar residual`           <dbl> 20.70, 1.60, 6.90, 8.50, 8.50, 6.90, 7.00,…
## $ cloruros                    <dbl> 0.045, 0.049, 0.050, 0.058, 0.058, 0.050, …
## $ `anhídrido sulfuroso libre` <dbl> 45, 14, 30, 47, 47, 30, 30, 45, 14, 28, 11…
## $ `anhídrido sulfuroso total` <dbl> 170, 132, 97, 186, 186, 97, 136, 170, 132,…
## $ densidad                    <dbl> 1.0010, 0.9940, 0.9951, 0.9956, 0.9956, 0.…
## $ pH                          <dbl> 3.00, 3.30, 3.26, 3.19, 3.19, 3.26, 3.18, …
## $ sulfatos                    <dbl> 0.45, 0.49, 0.44, 0.40, 0.40, 0.44, 0.47, …
## $ alcohol                     <dbl> 8.8, 9.5, 10.1, 9.9, 9.9, 10.1, 9.6, 8.8, …
## $ calidad                     <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 7, …
## $ variedad                    <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
set.seed(661) #fijo la semilla 

sample_var_1 = datos_tp %>% filter(datos_tp$variedad == 1)
sample_var_2 = datos_tp %>% filter(datos_tp$variedad == 2)
index_sample_data_1 <- sample(1:nrow(sample_var_1), size=1000, replace=FALSE)
index_sample_data_2 <- sample(1:nrow(sample_var_2), size=1000, replace=FALSE)
sample_data_var_1 = sample_var_1[index_sample_data_1,]
sample_data_var_2 = sample_var_2[index_sample_data_2,]
sample_data <- rbindlist(list(sample_data_var_1, sample_data_var_2))    # Rbind data.tables

glimpse(sample_data)
## Rows: 2,000
## Columns: 13
## $ `acidez fija`               <dbl> 5.9, 5.8, 9.1, 6.9, 7.2, 7.2, 7.4, 5.2, 7.…
## $ `acidez volátil`            <dbl> 0.320, 0.300, 0.280, 0.250, 0.230, 0.230, …
## $ `ácido cítrico`             <dbl> 0.26, 0.27, 0.49, 0.35, 0.32, 0.19, 0.27, …
## $ `azúcar residual`           <dbl> 1.5, 1.7, 2.0, 9.2, 8.5, 13.7, 1.2, 1.1, 1…
## $ cloruros                    <dbl> 0.057, 0.014, 0.059, 0.034, 0.058, 0.052, …
## $ `anhídrido sulfuroso libre` <dbl> 17.0, 45.0, 10.0, 42.0, 47.0, 47.0, 27.0, …
## $ `anhídrido sulfuroso total` <dbl> 141, 104, 112, 150, 186, 197, 99, 69, 112,…
## $ densidad                    <dbl> 0.99170, 0.98914, 0.99580, 0.99470, 0.9956…
## $ pH                          <dbl> 3.24, 3.40, 3.15, 3.21, 3.19, 3.12, 3.19, …
## $ sulfatos                    <dbl> 0.36, 0.56, 0.46, 0.36, 0.40, 0.53, 0.33, …
## $ alcohol                     <dbl> 10.700000, 12.600000, 10.100000, 11.500000…
## $ calidad                     <dbl> 5, 7, 5, 6, 6, 5, 6, 6, 5, 5, 7, 6, 6, 6, …
## $ variedad                    <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
print(paste("Balanced dataset - Variedad 1: ", nrow((sample_data %>% filter(variedad == 1)))))
## [1] "Balanced dataset - Variedad 1:  1000"
print(paste("Balanced dataset - Variedad 2: ", nrow((sample_data %>% filter(variedad == 2)))))
## [1] "Balanced dataset - Variedad 2:  1000"
Analizo el dataset
colnames(sample_data) <- c('acidez_fija','acidez_volatil','acido_citrico','azucar_residual','cloruros','anhidrido_sulfuroso_libre','anhidrido_sulfuroso_total','densidad','pH','sulfatos','alcohol','calidad','variedad')
sample_data
##       acidez_fija acidez_volatil acido_citrico azucar_residual cloruros
##    1:         5.9           0.32          0.26             1.5    0.057
##    2:         5.8           0.30          0.27             1.7    0.014
##    3:         9.1           0.28          0.49             2.0    0.059
##    4:         6.9           0.25          0.35             9.2    0.034
##    5:         7.2           0.23          0.32             8.5    0.058
##   ---                                                                  
## 1996:         5.6           0.54          0.04             1.7    0.049
## 1997:         7.3           0.65          0.00             1.2    0.065
## 1998:         6.9           0.41          0.31             2.0    0.079
## 1999:         9.0           0.38          0.41             2.4    0.103
## 2000:         9.2           0.67          0.10             3.0    0.091
##       anhidrido_sulfuroso_libre anhidrido_sulfuroso_total densidad   pH
##    1:                        17                       141  0.99170 3.24
##    2:                        45                       104  0.98914 3.40
##    3:                        10                       112  0.99580 3.15
##    4:                        42                       150  0.99470 3.21
##    5:                        47                       186  0.99560 3.19
##   ---                                                                  
## 1996:                         5                        13  0.99420 3.72
## 1997:                        15                        21  0.99460 3.39
## 1998:                        21                        51  0.99668 3.47
## 1999:                         6                        10  0.99604 3.13
## 2000:                        12                        48  0.99888 3.31
##       sulfatos alcohol calidad variedad
##    1:     0.36    10.7       5        1
##    2:     0.56    12.6       7        1
##    3:     0.46    10.1       5        1
##    4:     0.36    11.5       6        1
##    5:     0.40     9.9       6        1
##   ---                                  
## 1996:     0.58    11.4       5        2
## 1997:     0.47    10.0       7        2
## 1998:     0.55     9.5       6        2
## 1999:     0.58    11.9       7        2
## 2000:     0.54     9.5       6        2

1.- Aplique el Análisis de Componentes Principales a la base de datos. Presente los resultados y gráficos que considere adecuados. Interprete los resultados.

Me quedo con las variables numéricas para avanzar con el análisis
df_pca = sample_data
df_numericas_pca <-df_pca %>% dplyr::select(where(is.numeric))
Realizo un análisis univariado de cada atributo del conjunto de datos
data_long <- melt(df_pca)  
## Using variedad as id variables
#data_long
ggplot(data_long, aes(x=variable, y=value)) + 
    geom_boxplot() +
    facet_wrap(~variable, scale="free")

Realizo un análisis univariado de cada atributo del conjunto de datos, segregado por variedad de vinos
#Boxplot diferenciado por variedad de vino
ggplot(data_long, aes(x=variable, y=value, fill= variedad)) + 
    geom_boxplot() +
    facet_wrap(~variable, scale="free")

Genero una matriz de correlación entre las variables para el correspondiente análisis
#Matriz de correlación
m_cor <- cor(df_numericas_pca) 

# representa la matriz de correlaciones mediante círculos
corrplot(m_cor,method="circle")

Representacion de la varianza proporcionada por cada variable
variance_pca <- df_pca %>% 
  summarise_if(is.numeric, var) %>% 
  t() %>% 
  data.frame() %>% 
  round(3)

variance_pca
##                                  .
## acidez_fija                  2.436
## acidez_volatil               0.036
## acido_citrico                0.028
## azucar_residual             17.975
## cloruros                     0.002
## anhidrido_sulfuroso_libre  339.797
## anhidrido_sulfuroso_total 3715.930
## densidad                     0.000
## pH                           0.027
## sulfatos                     0.030
## alcohol                      1.352
## calidad                      0.761
Generación de componentes PCA para el conjunto de datos
pca <- prcomp(df_numericas_pca,scale = TRUE)# con datos estandarizados
names(pca)
## [1] "sdev"     "rotation" "center"   "scale"    "x"
Cargas (Loadings) para las componentes generadas
round(pca$rotation,2)
##                             PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9
## acidez_fija                0.27 -0.31  0.38 -0.25  0.14 -0.17  0.35 -0.20  0.20
## acidez_volatil             0.39  0.02 -0.31 -0.03  0.06  0.20  0.61  0.10  0.10
## acido_citrico             -0.13 -0.27  0.51 -0.02 -0.15 -0.40 -0.06  0.38  0.13
## azucar_residual           -0.31 -0.34 -0.14  0.09  0.47  0.19  0.03  0.43 -0.38
## cloruros                   0.31 -0.24  0.02  0.36 -0.42  0.40 -0.14  0.45  0.27
## anhidrido_sulfuroso_libre -0.40 -0.13 -0.13  0.38 -0.16 -0.13  0.42 -0.18  0.33
## anhidrido_sulfuroso_total -0.46 -0.15 -0.12  0.19 -0.16 -0.04  0.22 -0.05  0.01
## densidad                   0.22 -0.51 -0.11  0.07  0.42 -0.13  0.00  0.00  0.06
## pH                         0.20  0.30 -0.31  0.35  0.26 -0.59 -0.17  0.22  0.27
## sulfatos                   0.29 -0.10  0.22  0.63 -0.08 -0.10  0.03 -0.35 -0.55
## alcohol                   -0.01  0.46  0.33  0.06  0.06 -0.06  0.47  0.43 -0.25
## calidad                   -0.12  0.22  0.41  0.31  0.50  0.43 -0.09 -0.16  0.41
##                            PC10  PC11  PC12
## acidez_fija               -0.29 -0.32 -0.43
## acidez_volatil             0.52  0.17 -0.07
## acido_citrico              0.47  0.27  0.02
## azucar_residual           -0.09  0.09 -0.41
## cloruros                  -0.22 -0.16 -0.04
## anhidrido_sulfuroso_libre -0.34  0.44 -0.01
## anhidrido_sulfuroso_total  0.33 -0.72  0.08
## densidad                  -0.10 -0.05  0.68
## pH                        -0.01 -0.16 -0.24
## sulfatos                   0.11  0.05 -0.09
## alcohol                   -0.29 -0.13  0.32
## calidad                    0.20 -0.02  0.00
Correlación entre componentes generadas y atributos del conjunto de daots
contrib <- as.matrix(round(pca$rotation,2))
corrplot(contrib,is.corr=FALSE)

Represento la varianza explicada por cada componente y de forma acumulada
round(pca$center,2) #vector de medias
##               acidez_fija            acidez_volatil             acido_citrico 
##                      7.60                      0.40                      0.31 
##           azucar_residual                  cloruros anhidrido_sulfuroso_libre 
##                      4.52                      0.07                     25.93 
## anhidrido_sulfuroso_total                  densidad                        pH 
##                     93.24                      1.00                      3.25 
##                  sulfatos                   alcohol                   calidad 
##                      0.58                     10.47                      5.76
prop_varianza <- pca$sdev^2 / sum(pca$sdev^2)
prop_varianza #varianza explicada por cada componente
##  [1] 0.27683825 0.20521757 0.15634119 0.08242453 0.06948051 0.04950167
##  [7] 0.04356408 0.03971796 0.03753468 0.02186397 0.01427615 0.00323943
prop_varianza_acum_pca <- cumsum(prop_varianza)
prop_varianza_acum_pca #varianza acumilada explicada por los componentes de forma creciente segun su importancia
##  [1] 0.2768383 0.4820558 0.6383970 0.7208215 0.7903021 0.8398037 0.8833678
##  [8] 0.9230858 0.9606205 0.9824844 0.9967606 1.0000000
Representación gráfica de la varianza explicada por cada componente y de forma acumulada
ggplot(data = data.frame(prop_varianza_acum_pca, pc = 1:12),
       aes(x = pc, y = prop_varianza_acum_pca, group = 1)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  labs(x = "Componente principal",
       y = "Prop. varianza explicada acumulada")

screeplot(pca, type = "l", npcs = 12)
abline(h = 1, col="red", lty=5)
legend("topright", legend=c("Eigenvalue = 1"),
       col=c("red"), lty=5, cex=0.6)

fviz_eig(pca, ncp =12, addlabels = TRUE, main="")

Visualización de scores para todas las componentes principales
#pca_2$x[,1]# scores para la primer componente (PC1)
res.ind <- get_pca_ind(pca)
head(res.ind$coord,10)
##          Dim.1      Dim.2       Dim.3      Dim.4        Dim.5       Dim.6
## 1  -1.01189442  1.3191155 -0.77274712 -1.0357319 -1.303163205 -0.04051432
## 2  -1.66875699  3.3360841  0.67248805  0.9158157  0.003178629 -0.32182485
## 3   0.01602419 -0.8192638  0.82974834 -1.5371725 -0.683656840 -0.75875675
## 4  -2.33345523  0.2698171  0.04272720 -0.2733312  0.534078552 -0.07713830
## 5  -2.33834010 -0.7957860 -0.42077524  0.1231915  0.134039028  0.15845606
## 6  -2.23501071 -2.3477902 -1.56768877  0.2210427  0.509138131  0.22625340
## 7  -1.14291956  0.7835540  0.05342973 -1.2018491 -0.456717256  0.13685204
## 8  -0.71628256  1.6486190 -0.14371084 -0.3971398 -0.535746151  0.10267295
## 9  -0.51483343  0.4081233 -0.23681707 -1.5712144 -1.066796279  0.21658374
## 10 -0.68431517  1.1294388  0.37372431 -0.7401340 -0.808225877 -1.52074829
##           Dim.7       Dim.8       Dim.9      Dim.10     Dim.11       Dim.12
## 1  -0.495919296  0.35113630 -0.05141232  0.25368656 -0.6094454  0.150398260
## 2   0.470909358 -0.13683111  0.23585111 -0.06141238  0.3736711 -0.280061337
## 3  -0.391449922  0.06641265  0.11350261  0.21009852 -0.6343479  0.119948769
## 4   0.435569228  0.76019151  0.04545350 -0.41840389 -0.1270732  0.180397465
## 5  -0.009342246  0.04604676  0.55799854 -0.22347452 -0.4765575 -0.015741529
## 6  -0.035025845 -0.30334927 -0.77677646 -0.63039045 -0.5180117  0.009249802
## 7  -0.815847371 -0.65442870  0.89109712 -0.22575357 -0.1096548 -0.145302705
## 8  -1.296704575 -0.60975149 -0.23873680  0.62686972  0.6294417  0.047023266
## 9  -0.786231545 -0.47653397 -0.34558505 -0.01493300 -0.5615664 -0.048033200
## 10 -1.155540815  0.44736925 -0.28989256  0.03669330 -0.1287600 -0.015182979
Grafico biplot para comprender la explicación de cada componente generada según la variedad y los autovectores de cada atributo
autoplot(pca, 
         data = df_pca, 
         colour = 'variedad',
         loadings = TRUE, 
         loadings.colour = 'black',
         loadings.label = TRUE, 
         loadings.label.size = 3)

Conclusiones

Podemos obtener como resultado del PCA, que sería prudente quedarnos con las primeras 7 componentes que representan más del 86% de la variabilidad de la información, además se comienza a aplanar la curva que explica la diferencia de vairabilidad entre las siguientes componentes

Observamos que los principales atributos con los que nos quedamos a posteriri de decidir con la información de las componentes serian anhidrido_sulfuroso_total (1) densidad (2) acido_citrico (3) sulfatos (4) calidad (5) pH (6) alcohol (7)

anhidrido_sulfuroso_total y anhidrido_sulfuroso_libre correlacionan de forma casi total

Mientras que otras variables tienne una correlación alta dada por el ángulo que forman sus autovectores Azucar_resicual/Acido_citrico, acidez_fija/cloruros, sulfatos/acidez_volatil

Por otro lado, vemos que otras variables no tienen ningún tipo de correlación y las podemos determinar independientes entre si, como alcohol/acidez_volatil, ph/calidad, acidez_fija/azucar_residual

2.- Realice el Análisis Discriminante para clasificar los vinos según la variable variedad de vino. Interprete los resultados.

#Separo el dataset para trabajar el analisis discriminante
df_ad = sample_data
Grafico la distribución de cada variable separada por variedad para su correspondiente análisis
ggpairs(df_ad, legend = 1, columns = 1:12, aes(color = variedad, alpha = 0.5),
        upper = list(continuous = "points"))+
  theme(legend.position = "bottom")

Separo el dataset en entrenamiento y prueba para poder realizar la predicción según el análisis discriminate y no sobreajustar los resultados de la predicción
#Separo en train y test para predecir

set.seed(661)#setear la semilla
# Create data split for train and test
df_split_ad <- initial_split(df_ad,
                          prop = 0.9,
                          strata = variedad)#en este caso no es necesario porque las clases están balanceadas

df_train_ad <- df_split_ad %>% training()

df_test_ad <- df_split_ad %>% testing()

# Número de datos en test y train
paste0("Total del dataset de entrenamiento: ", nrow(df_train_ad))
## [1] "Total del dataset de entrenamiento: 1800"
paste0("Total del dataset de testeo: ", nrow(df_test_ad))
## [1] "Total del dataset de testeo: 200"
Genero sub conjuntos por variedad de vinos para hacer el correspondiente análisis de normalidad y homocedasticidad
blanco_ad<- subset(df_ad[,1:12], df_ad$variedad == "1")
tinto_ad <- subset(df_ad[,1:12], df_ad$variedad == "2")
Valido normalidad con un test de Shapiro-Wilk
mvShapiro.Test(as.matrix(blanco_ad))
## 
##  Generalized Shapiro-Wilk test for Multivariate Normality by
##  Villasenor-Alva and Gonzalez-Estrada
## 
## data:  as.matrix(blanco_ad)
## MVW = 0.91982, p-value < 2.2e-16
Valido normalidad con un test de Shapiro-Wilk
#Valido normalidad => no se cumple porque es < 0,05
mvShapiro.Test(as.matrix(tinto_ad))
## 
##  Generalized Shapiro-Wilk test for Multivariate Normality by
##  Villasenor-Alva and Gonzalez-Estrada
## 
## data:  as.matrix(tinto_ad)
## MVW = 0.88076, p-value < 2.2e-16
Valido homocedasticidad con boxM
#Matriz vairanza y covarianza son iguales, valido homocedasticidad => no se cumple tampoco

boxM(data = df_ad[, 1:12], grouping = df_ad$variedad)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  df_ad[, 1:12]
## Chi-Sq (approx.) = 4289.5, df = 78, p-value < 2.2e-16
#Por mas que nada se cumple seguimos el analisis, sabiendo que los resultados pueden no tener la validez que esperamos

Conclusiones

Determino independencia suponiendo que los datos se obtuvieron de fuentes independientes

Vemos que no se cumple la normalidad en ninguna de las dos distribuciones planteadas

Vemos que no se cumple homocedasticidad en el conjunto de datos

Por mas que los supuestos no se cumplen, avanzamos con el análisis a sabiendas que los resultados pueden no tener la validez que esperamos

Ejecuto el análisis de discriminantes para obtener la función de dimensiones que nos separen los conjuntos de datos
model_lda <- lda(variedad~., data=df_train_ad)
model_lda
## Call:
## lda(variedad ~ ., data = df_train_ad)
## 
## Prior probabilities of groups:
##   1   2 
## 0.5 0.5 
## 
## Group means:
##   acidez_fija acidez_volatil acido_citrico azucar_residual   cloruros
## 1    6.860278      0.2821444     0.3411222        6.524278 0.04621778
## 2    8.372889      0.5232333     0.2759889        2.574833 0.08808333
##   anhidrido_sulfuroso_libre anhidrido_sulfuroso_total  densidad       pH
## 1                  36.08000                 140.10167 0.9941380 3.186922
## 2                  15.94556                  46.83222 0.9967918 3.305744
##    sulfatos  alcohol  calidad
## 1 0.4923333 10.47361 5.858889
## 2 0.6612444 10.43659 5.646667
## 
## Coefficients of linear discriminants:
##                                    LD1
## acidez_fija                -0.44594457
## acidez_volatil              2.05824631
## acido_citrico              -0.58935312
## azucar_residual            -0.36366284
## cloruros                    2.40213219
## anhidrido_sulfuroso_libre   0.01517737
## anhidrido_sulfuroso_total  -0.01876167
## densidad                  898.41125911
## pH                         -1.19965511
## sulfatos                    0.92102307
## alcohol                     0.75018119
## calidad                     0.05713148
Represento la varianza entre grupos explicada por cada FD, en este caso hay solo una que representa prácticamente el 100% de la varianza
prop = model_lda$svd^2/sum(model_lda$svd^2)
prop 
## [1] 1
#Coeficiente lda calculado para cada entrada del set analizado
lda.data <- cbind(df_train_ad, predict(model_lda)$x)
#lda.data
#ggplot(lda.data, aes(LD1)) +
#  geom_point(aes(color = variedad))
#Predicciones sobre el set de entrenamiento para cada entrada del set analizado
predictions_train <- model_lda %>% predict(df_train_ad)
#predictions_train$class
Matriz de confusión con los resultados de la predicción sobre el conjunto de entrenamiento
table(predict(model_lda,type="class")$class,df_train_ad$variedad)
##    
##       1   2
##   1 898  12
##   2   2 888
#partimat (variedad~. , data=df_train_ad , method="lda")
#ggord(model_lda, df_train_ad$variedad)
predictions <- model_lda %>% predict(df_test_ad)
predictions$class
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [186] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Levels: 1 2
Matriz de confusión con los resultados de la predicción sobre el conjunto de prueba
lda.test <- predict(model_lda,df_test_ad)
df_test_ad$lda <- lda.test$class
table(df_test_ad$lda,df_test_ad$variedad)
##    
##       1   2
##   1 100   0
##   2   0 100

Si los supuesto no se cumplen, podrían no ser válidos los resultados obtenidos Como se utiliza las matrices de varianzas y covarianzas, estas técnicas son sensibles a outliers

model_qda <- qda(variedad ~ ., df_train_ad)
model_qda
## Call:
## qda(variedad ~ ., data = df_train_ad)
## 
## Prior probabilities of groups:
##   1   2 
## 0.5 0.5 
## 
## Group means:
##   acidez_fija acidez_volatil acido_citrico azucar_residual   cloruros
## 1    6.860278      0.2821444     0.3411222        6.524278 0.04621778
## 2    8.372889      0.5232333     0.2759889        2.574833 0.08808333
##   anhidrido_sulfuroso_libre anhidrido_sulfuroso_total  densidad       pH
## 1                  36.08000                 140.10167 0.9941380 3.186922
## 2                  15.94556                  46.83222 0.9967918 3.305744
##    sulfatos  alcohol  calidad
## 1 0.4923333 10.47361 5.858889
## 2 0.6612444 10.43659 5.646667
Matriz de confusión con los resultados de la predicción sobre el conjunto de entrenamiento evaluado con QDA
table(predict(model_qda,type="class")$class,df_train_ad$variedad)
##    
##       1   2
##   1 882  12
##   2  18 888
Matriz de confusión con los resultados de la predicción sobre el conjunto de prueba evaluado con QDA
lda.test_qda <- predict(model_qda,df_test_ad)
df_test_ad$qda <- lda.test_qda$class
table(df_test_ad$qda,df_test_ad$variedad)
##    
##       1   2
##   1 100   0
##   2   0 100

Conclusiones Podemos determinar a partir de la predicción y el analisis de su matriz de confusión, la separación generada por LDA en el conjunto de datos, arroja un accuracy del 100% en el dataset de prueba

Solo tenemos una dimensión debido a que nuestras variables objetivos solamente tienen dos posibilidades

Por el tipo de dato a validar y el dataset analizado, podría tener mayor sentido analizar otra variable para separa el mismo como por ejemplo la calidad donde nos arrojarían varias dimensiones a separar

En principio para el objetivo del análisis, determinamos que LDA funciona muy bien en su predicción de vinos tintos y blancos

No vemos una mejora significativa utilizando QDA

3.- Aplique el algoritmo SVM al conjunto de datos. Interprete los resultados.

df_svm = sample_data
Analizo correlación de datos
#Veo como correlacionan los datos
df_svm %>% ggpairs(.,mapping=ggplot2::aes(color = variedad,alpha = 0.1),
        upper = list(continuous = wrap("cor", size = 2.5),discrete = "blank", combo="blank"),
        lower = list(combo = "box"),progress = F)+
        theme+
        labs(title= 'Descripción de variables en la base de datos', x='Variable', y='Variable')+
        scale_fill_manual(values=c('royalblue2','red'))+
        scale_color_manual(values=c('royalblue2','#ff7474ff'))

# Grafico biplot de pca para poder comparar posteriormente con analisis de svm
datos_para_acp = df_svm %>% 
      dplyr::select(is.numeric) # todas las variables numéricas
## Warning: Predicate functions must be wrapped in `where()`.
## 
##   # Bad
##   data %>% select(is.numeric)
## 
##   # Good
##   data %>% select(where(is.numeric))
## 
## ℹ Please update your code.
## This message is displayed once per session.
df_svm.pc = prcomp(datos_para_acp,scale = TRUE) #escalo los datos
#grafico
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=0.5,
         groups=factor(df_svm$variedad)) +
        scale_color_manual(name="variedad",values=c('royalblue2','#ff7474ff'),
                           labels=c("blanco",'tinto')) +
        theme + labs(title='Análisis de componentes principales') + 
        theme(legend.position=c(.85,.15)) 

#Analisis estadistico univariado agrupado por variedad => Estaria bueno tenerlo pero no lo puedo hacer andar de momento

# Creo dataframes de acuerdo al diagnóstico
 data_n <- df_svm%>%filter(variedad==1)
 data_m <- df_svm%>%filter(variedad==2)
Analizo normalidad multivariada y multivariada por grupo con Shapiro Wilks
pval_all <- mshapiro.test(t(df_svm[,1:12]))
pval_blanco <- mshapiro.test(t(df_svm%>%dplyr::filter(variedad==1)%>% dplyr::select(1:12)))
pval_tinto <- mshapiro.test(t(df_svm%>%dplyr::filter(variedad==2)%>% dplyr::select(1:12)))

ShapiroW_p.valor <- c(pval_all$p.value, pval_blanco$p.value, pval_tinto$p.value)
Datos <- c('todos','subset blanco', 'subset tinto')

# Imprimo resultados de normalidad multivariada total y por grupos (estos últimos son relevantes)
kable(cbind(Datos,ShapiroW_p.valor))
Datos ShapiroW_p.valor
todos 5.24820022494157e-43
subset blanco 1.32661116153684e-33
subset tinto 6.99150814871208e-43
Analizo igualdad de matrices de varianzas y covarianzas con boxM
p.valor <- boxM(data = df_svm[, 1:12], grouping = df_svm$variedad)$p.value
Test <- boxM(data = df_svm[, 1:12], grouping = df_svm$variedad)$method

# Imprimo resultados de test
kable(cbind(Test, p.valor))
Test p.valor
Box’s M-test for Homogeneity of Covariance Matrices 0
Como boxM es sensible a la falta de normalidad, aplico Levene utilizando betadisper del paquete “vegan” (equivalente a levene, pero multivariado)
matriz_de_distancias <- vegan::betadisper(dist(df_svm[, 1:12], method='euclidean'), df_svm$variedad, type = c("median","centroid"), bias.adjust = T,sqrt.dist = FALSE, add = FALSE)
## Registered S3 methods overwritten by 'vegan':
##   method      from      
##   rev.hclust  dendextend
##   plot.rda    klaR      
##   predict.rda klaR      
##   print.rda   klaR
test_levene <- anova(matriz_de_distancias)
p.valor <- test_levene$`Pr(>F)`[1]
TukeyHSD(matriz_de_distancias)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = distances ~ group, data = df)
## 
## $group
##          diff       lwr       upr p adj
## 2-1 -12.98251 -15.28052 -10.68451     0
 #plot(matriz_de_distancias)
Test <- 'Levene'

# Imprimo resultado del test
kable(cbind(Test,p.valor))
Test p.valor
Levene 9.94794504885316e-28
Analizo cómo me da Hotelling para ver diferencias en el vector de medias de cada grupo
HOTELLING <- HotellingsT2Test(as.matrix(df_svm[,-c(13)]) ~ variedad, data =df_svm)
Test <- HOTELLING$method
p.valor <- HOTELLING$p.value[1]

# Imprimo resultados en una tabla
kable(cbind(Test, p.valor ))
Test p.valor
Hotelling’s two sample T2-test 0 
# se utiliza el paquete npmv no paramétrico para comparar vector de medias. (Nonparametric Inference for Multivariate Data: R Package npmv, January 2017, Volume 76, Issue 4. doi: 10.18637/jss.v076.i04, https://www.jstatsoft.org/article/view/v076i04)
# => Agregar todas las variables
#noparam <- nonpartest(pH | calidad | densidad | alcohol | sulfatos | cloruros ~ variedad, data = df_svm, permreps = 1000, plots=F) 
noparam <- nonpartest(  acidez_fija | acidez_volatil | acido_citrico | azucar_residual | cloruros | anhidrido_sulfuroso_libre | anhidrido_sulfuroso_total | densidad | pH | sulfatos | alcohol | calidad ~ variedad, data = df_svm, permreps = 1000, plots=F) 
p.valor <- noparam$results$`P-value`[1]
Test <- 'No paramétrico multivariado (npmv)'

# Imprimo el resultado
kable(cbind(Test, p.valor ))
Test p.valor
No paramétrico multivariado (npmv) 0
Separo conjunto de entrenamiento y prueba para una posterior predicción evitando resultado sobreajustados solo sobre el conjunto de entrenamiento
set.seed(661) # para asegurar reproducibilidad
dt = sort(sample(nrow(df_svm), nrow(df_svm)*.8))
datos_train_svm<-df_svm[dt,]
datos_test_svm<-df_svm[-dt,]
Se realiza el escalado/estandarización con -> (x - mean(x)) / sd(x)
# Se realiza el escalado/estandarización con ->  (x - mean(x)) / sd(x)

# Calculo media y sd de subconjunto de entrenamiento (train), y con esos datos hago el escalado del test. La idea de escalar el conjunto de test (prueba) utilizando datos solamente de train es para evitar el data leakeage.

# Hago escalado a mano del test set con media del training, y sd del training
for (k in 1:12){
  datos_test_svm[,k]=round(print(datos_test_svm[,..k]-colMeans(datos_test_svm[,..k]))/sd(unlist(datos_train_svm[,..k])),2)
}
##      acidez_fija
##   1:    -0.53625
##   2:    -0.33625
##   3:    -1.13625
##   4:    -1.13625
##   5:    -2.73625
##  ---            
## 396:    -0.43625
## 397:    -0.03625
## 398:     3.16375
## 399:    -1.63625
## 400:    -0.83625
##      acidez_volatil
##   1:     -0.1621625
##   2:     -0.1821625
##   3:      0.2278375
##   4:     -0.0521625
##   5:     -0.0421625
##  ---               
## 396:     -0.0721625
## 397:      0.3128375
## 398:     -0.1821625
## 399:      0.1678375
## 400:      0.0178375
##      acido_citrico
##   1:      -0.12185
##   2:      -0.04185
##   3:      -0.11185
##   4:       0.02815
##   5:      -0.06185
##  ---              
## 396:      -0.08185
## 397:      -0.21185
## 398:       0.17815
## 399:      -0.31185
## 400:      -0.00185
##      azucar_residual
##   1:        9.123625
##   2:       -3.376375
##   3:        4.123625
##   4:       -1.976375
##   5:        3.223625
##  ---                
## 396:       -2.276375
## 397:       -1.976375
## 398:       -1.776375
## 399:       -2.376375
## 400:       -2.576375
##       cloruros
##   1: -0.014905
##   2: -0.025905
##   3: -0.020905
##   4: -0.015905
##   5: -0.035905
##  ---          
## 396: -0.000905
## 397:  0.017095
## 398:  0.021095
## 399:  0.012095
## 400:  0.012095
##      anhidrido_sulfuroso_libre
##   1:                    21.765
##   2:                     1.765
##   3:                    55.765
##   4:                    15.265
##   5:                    -1.235
##  ---                          
## 396:                     9.765
## 397:                   -16.235
## 398:                   -14.235
## 399:                   -19.235
## 400:                    -4.235
##      anhidrido_sulfuroso_total
##   1:                 102.85625
##   2:                   4.85625
##   3:                 129.85625
##   4:                 115.85625
##   5:                  21.85625
##  ---                          
## 396:                 -24.14375
## 397:                 -68.14375
## 398:                 -62.14375
## 399:                 -85.14375
## 400:                 -43.14375
##         densidad
##   1:  0.00306995
##   2: -0.00288005
##   3:  0.00046995
##   4: -0.00244005
##   5: -0.00317005
##  ---            
## 396:  0.00029995
## 397:  0.00201995
## 398:  0.00161995
## 399: -0.00078005
## 400:  0.00109995
##           pH
##   1: -0.1177
##   2: -0.0477
##   3: -0.0677
##   4: -0.0877
##   5:  0.1523
##  ---        
## 396:  0.1923
## 397:  0.1523
## 398: -0.0177
## 399:  0.3523
## 400:  0.2323
##       sulfatos
##   1: -0.033625
##   2: -0.233625
##   3: -0.123625
##   4:  0.046375
##   5: -0.163625
##  ---          
## 396:  0.056375
## 397: -0.073625
## 398:  0.116375
## 399: -0.023625
## 400: -0.013625
##         alcohol
##   1: -1.4245417
##   2: -0.6245417
##   3: -1.1245417
##   4: -0.4245417
##   5:  0.8754583
##  ---           
## 396: -0.3245417
## 397: -0.7245417
## 398:  1.2754583
## 399:  1.0754583
## 400: -0.9245417
##      calidad
##   1: -0.7375
##   2:  0.2625
##   3: -0.7375
##   4: -0.7375
##   5:  0.2625
##  ---        
## 396: -0.7375
## 397: -0.7375
## 398:  0.2625
## 399:  0.2625
## 400:  0.2625
# Hago automáticamente con la función scale, el escalado de training
datos_train_svm = as.data.frame(scale(datos_train_svm[,1:12]))
datos_train_svm$variedad <- df_svm[dt,]$variedad
# y las pongo en = orden que testset
# creo una única df con todos los datos escalados (que usaré luego para clustering)
datos_escalados <- as.data.frame(scale(df_svm[,1:12]))
datos_escalados$variedad <- df_svm$variedad
Defino modelo SVM lineal, separo los conjuntos de datos para entrenamiento y prueba, recolecto las métricas que me ayudaran a sacar conclusiones sobre las predicciones y grafico la curva ROC
# Defino modelo SVM => Tengo que agregar todas las columnas
sample_svn_tr = datos_train_svm[,1:13]
sample_svn_te = datos_test_svm[,1:13]

set.seed(661)
task = makeClassifTask(data = sample_svn_tr, target = "variedad") 
lrn_svm1 = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost=2)) 
mod_svm1 = mlr::train(lrn_svm1, task)

# Predicción TEST
pred_svm1= predict(mod_svm1, newdata = sample_svn_te)
## Warning in predict.WrappedModel(mod_svm1, newdata = sample_svn_te): Provided
## data for prediction is not a pure data.frame but from class data.table, hence it
## will be converted.
acc_svm1 <- round(measureACC(pred_svm1$data$truth, pred_svm1$data$response),3)
AUC_svm1_te <- round(measureAUC(as.data.frame(pred_svm1)$prob.2,as.data.frame(pred_svm1)$truth,1,2),3)

# ···························································································
# Predicción TRAIN (naive)
pred_svm_1 = predict(mod_svm1, newdata = sample_svn_tr) # para comparar con set de test
acc_svm_1 <- round(measureACC(as.data.frame(pred_svm_1)$truth, as.data.frame(pred_svm_1)$response),3)
AUC_svm1_tr <- round(measureAUC(as.data.frame(pred_svm_1)$prob.2,as.data.frame(pred_svm_1)$truth, 1,2),3)

#print(sample_svn_tr)
# ···························································································
# Cambio el threshold [esto lo hago para train y test]
acc=NULL
acc2=NULL
threshold = seq(0.1,0.95,0.01)
for (i in 1:length(threshold)) {
        pred = setThreshold(pred_svm1, threshold = threshold[i])
        acc[i] = measureACC(as.data.frame(pred)$truth, as.data.frame(pred)$response)}
for (i in 1:length(threshold)) {
        pred2 = setThreshold(pred_svm_1, threshold = threshold[i])
        acc2[i] = measureACC(as.data.frame(pred2)$truth, as.data.frame(pred2)$response)}
par(mfrow = c(2,5))

new_df1 <- as.data.frame(cbind(threshold,acc))
new_df1 <- new_df1%>%mutate(sub_data='test')
new_df2 <- as.data.frame(cbind(threshold,acc2))
colnames(new_df2) <- c('threshold','acc')
new_df2 <- new_df2%>%mutate(sub_data='train')

new_df <- as.data.frame(rbind(new_df1,new_df2))

#new_df1[which.max(new_df1$acc),"threshold"] 
#new_df2[which.max(new_df2$acc),"threshold"] 
# ···························································································
# Gráfico de cómo varía la métrica de performance accuracy, de acuerdo al umbral elegido
ggplot(new_df, aes(x=threshold, y=acc)) + geom_line(aes(color = sub_data,linetype=sub_data)) +
        theme +labs(x='Umbral', y='Métrica de performance (accuracy)', 
                     title= 'Evaluación del modelo de Máquinas de soporte vectorial SVM') +
        scale_color_manual(values = c("red", "darkred"),labels=c('prueba','entrenamiento')) +
        scale_linetype_manual(values=c(1,2), labels=c('prueba','entrenamiento')) + 
        labs(color='Conjunto de\n evaluación',linetype='Conjunto de\n evaluación')

#Representamos la matriz de confusion para poder hacer un analisis de las diferentes metricas sobre la prediccion

#install.packages("caret")
library(caret)
## 
## Attaching package: 'caret'
## The following objects are masked from 'package:DescTools':
## 
##     MAE, RMSE
## The following object is masked from 'package:mlr':
## 
##     train
## The following object is masked from 'package:purrr':
## 
##     lift
confusionMatrix(as.data.frame(pred_svm1)$truth, as.data.frame(pred_svm1)$response)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   1   2
##          1 199   0
##          2   2 199
##                                           
##                Accuracy : 0.995           
##                  95% CI : (0.9821, 0.9994)
##     No Information Rate : 0.5025          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.99            
##                                           
##  Mcnemar's Test P-Value : 0.4795          
##                                           
##             Sensitivity : 0.9900          
##             Specificity : 1.0000          
##          Pos Pred Value : 1.0000          
##          Neg Pred Value : 0.9900          
##              Prevalence : 0.5025          
##          Detection Rate : 0.4975          
##    Detection Prevalence : 0.4975          
##       Balanced Accuracy : 0.9950          
##                                           
##        'Positive' Class : 1               
##
Para independizarnos de la elección del umbral, grafico curvas ROC para las predicciones del modelo SVM con los datos de Entrenamiento y Prueba
df_svm_threshold = generateThreshVsPerfData(list(svm_te = pred_svm1, svm_tr = pred_svm_1), 
                                  measures = list(fpr, tpr, mmce))

plotROCCurves(df_svm_threshold) + theme +
        labs(title='Curva ROC del modelo de Máquinas de soporte vectorial SVM kernel lineal', 
             x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)',
             color='Conjunto de\n evaluación') +
        scale_color_manual(values = c("red", "darkred"), labels=c('prueba','entrenamiento')) +
        geom_label(label='AUC test', x=0.35, y=0.75, label.size = 0.3, size=4,
                   color = "red",fill="white") + 
        geom_label(label='AUC train',x=0.07, y=0.97, label.size = 0.3, size=4,
                   color = "darkred",fill="white")

# ················ Métricas del modelo de SVM ················
Métrica <- c('valor','datos')
Accuracy <- c(acc_svm1,'prueba')
Accuracy. <- c(acc_svm_1,'entrenamiento')
AUC_ROC <- c(AUC_svm1_te,'prueba')
AUC_ROC. <- c(AUC_svm1_tr,'entrenamiento')
# Imprimo resultados
kable(rbind(Métrica, Accuracy, Accuracy., AUC_ROC, AUC_ROC.))
Métrica valor datos
Accuracy 0.995 prueba
Accuracy. 0.992 entrenamiento
AUC_ROC 0.998 prueba
AUC_ROC. 0.996 entrenamiento
Buscando un modelo más robusto, defino modelo SVM radial, separo los conjuntos de datos para entrenamiento y prueba, recolecto las métricas que me ayudarán a sacar conclusiones sobre las predicciones y grafico la curva ROC
# Defino modelo SVM radial

sample_svn_tr_radial = datos_train_svm[,1:13]
sample_svn_te_radial = datos_test_svm[,1:13]

set.seed(661)
task = makeClassifTask(data = sample_svn_tr_radial, target = "variedad") 
lrn_svm3 = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "radial", cost=2)) 
mod_svm3 = mlr::train(lrn_svm3, task)
# Predicción TEST
pred_svm3= predict(mod_svm3, newdata = sample_svn_te_radial)
## Warning in predict.WrappedModel(mod_svm3, newdata = sample_svn_te_radial):
## Provided data for prediction is not a pure data.frame but from class data.table,
## hence it will be converted.
acc_svm3 <- round(measureACC(pred_svm3$data$truth, pred_svm3$data$response),3)
AUC_svm3_te <- round(measureAUC(as.data.frame(pred_svm3)$prob.2,as.data.frame(pred_svm3)$truth,1,2),3)
# ···························································································
# Predicción TRAIN (naive)
pred_svm_3 = predict(mod_svm3, newdata = sample_svn_tr_radial) # por si quiero ver naive sobre training
acc_svm_3 <- round(measureACC(as.data.frame(pred_svm_3)$truth, as.data.frame(pred_svm_3)$response),3)
AUC_svm3_tr <- round(measureAUC(as.data.frame(pred_svm_3)$prob.2,as.data.frame(pred_svm_3)$truth, 1,2),3)
# ···························································································
# Cambio el threshold [esto lo hago para train y test]
acc=NULL
acc2=NULL
threshold = seq(0.1,0.95,0.01)
for (i in 1:length(threshold)) {
        pred = setThreshold(pred_svm3, threshold = threshold[i])
        acc[i] = measureACC(as.data.frame(pred)$truth, as.data.frame(pred)$response)}
for (i in 1:length(threshold)) {
        pred2 = setThreshold(pred_svm_3, threshold = threshold[i])
        acc2[i] = measureACC(as.data.frame(pred2)$truth, as.data.frame(pred2)$response)}
par(mfrow = c(2,5))

new_df1 <- as.data.frame(cbind(threshold,acc))
new_df1 <- new_df1%>%mutate(sub_data='test')
new_df2 <- as.data.frame(cbind(threshold,acc2))
colnames(new_df2) <- c('threshold','acc')
new_df2 <- new_df2%>%mutate(sub_data='train')

new_df <- as.data.frame(rbind(new_df1,new_df2))

#new_df1[which.max(new_df1$acc),"threshold"] 
#new_df2[which.max(new_df2$acc),"threshold"] 
# ···························································································
# Gráfico de cómo varía la métrica de performance accuracy, de acuerdo al umbral elegido
ggplot(new_df, aes(x=threshold, y=acc)) + geom_line(aes(color = sub_data,linetype=sub_data)) +
        theme +labs(x='Umbral', y='Métrica de performance (accuracy)', 
                     title= 'Evaluación del modelo de Máquinas de soporte vectorial SVM') +
        scale_color_manual(values = c("red", "darkred"),labels=c('prueba','entrenamiento')) +
        scale_linetype_manual(values=c(1,2), labels=c('prueba','entrenamiento')) + 
        labs(color='Conjunto de\n evaluación',linetype='Conjunto de\n evaluación')

# Para independizarnos de la elección del umbral, grafico curvas ROC para las predicciones del modelo SVM con los datos de TEST y TRAIN
df_svm_threshold_radial = generateThreshVsPerfData(list(svm_te = pred_svm3, svm_tr = pred_svm_3), 
                                  measures = list(fpr, tpr, mmce))

plotROCCurves(df_svm_threshold_radial) 

confusionMatrix(as.data.frame(pred_svm3)$truth, as.data.frame(pred_svm3)$response)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   1   2
##          1 199   0
##          2   2 199
##                                           
##                Accuracy : 0.995           
##                  95% CI : (0.9821, 0.9994)
##     No Information Rate : 0.5025          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.99            
##                                           
##  Mcnemar's Test P-Value : 0.4795          
##                                           
##             Sensitivity : 0.9900          
##             Specificity : 1.0000          
##          Pos Pred Value : 1.0000          
##          Neg Pred Value : 0.9900          
##              Prevalence : 0.5025          
##          Detection Rate : 0.4975          
##    Detection Prevalence : 0.4975          
##       Balanced Accuracy : 0.9950          
##                                           
##        'Positive' Class : 1               
## 
# ················ Métricas del modelo de SVM ················
Métrica <- c('valor','datos')
Accuracy <- c(acc_svm3,'prueba')
Accuracy. <- c(acc_svm_3,'entrenamiento')
AUC_ROC <- c(AUC_svm3_te,'prueba')
AUC_ROC. <- c(AUC_svm3_tr,'entrenamiento')
# Imprimo resultados
kable(rbind(Métrica, Accuracy, Accuracy., AUC_ROC, AUC_ROC.))
Métrica valor datos
Accuracy 0.995 prueba
Accuracy. 0.995 entrenamiento
AUC_ROC 0.996 prueba
AUC_ROC. 0.999 entrenamiento
Realizo una comparación de las predicciones realizadas por ambos métodos SVM en contraste con la evaluación de componentes principales generadas anteriormente
SVM_metrics1 <- calculateROCMeasures(pred_svm1)
SVM_metrics3 <- calculateROCMeasures(pred_svm3)

# Ingenua para ver cómo le va 
pred_todos=NULL
pred_todos_svm1 <- as.data.frame(predict(mod_svm1, newdata = datos_escalados[-13]))
pred_todos_svm3 <- as.data.frame(predict(mod_svm3, newdata = datos_escalados[-13]))

# ······················ PCA datos originales ··················································
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=0.5,groups=factor(df_pca$variedad)) +
        scale_color_manual(name="variedad", values=c('royalblue2','#ff7474ff'),
                           labels=c("blanco","tinto")) +
        theme + labs(
        title= 'Representación de las etiquetas reales\nen las componentes principales 1 y 2')+
        theme(legend.position=c(.865,.15))

# ······················ PCA + predicciones SVM lineal (modelo 1) ·······························
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=.5,groups=factor(pred_todos_svm1$response)) +
        scale_color_manual(name="Predicción", values=c('royalblue2','#ff7474ff'),
                           labels=c("blanco","tinto")) +
        theme + labs(title= 'Representación de las predicciones ingenuas del modelo SVM (l)\nen las componentes principales 1 y 2') +
        theme(legend.position=c(.865,.15))

# ······················ PCA + predicciones SVM radial (modelo 3) ·······························
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=.5,groups=factor(pred_todos_svm3$response)) +
        scale_color_manual(name="Predicción", values=c('royalblue2','#ff7474ff'),
                           labels=c("blanco","tinto")) +
        theme + labs(title= 'Representación de las predicciones ingenuas del modelo SVM (r)\nen las componentes principales 1 y 2') +
        theme(legend.position=c(.865,.15))

# ······················ curvas ROC para todos los modelos ········································
df_todos = generateThreshVsPerfData(list(svm1 = pred_svm1,
                                         svm3 = pred_svm3), measures = list(fpr, tpr, mmce))

plotROCCurves(df_todos) + theme + 
        labs(title='Curvas ROC de modelos de clasificación supervisada (datos de prueba)',
             x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)', 
             color=' Modelo en\n evaluación') +
        scale_color_manual(values = c("red", "black"),
                           labels=c('SVM (l)','SVM (r)'))+
        theme(legend.position=c(0.915,0.25))

Conclusiones Podemos determinar a partir del análisis realizado, de que para este caso, SVM está realizando una separación muy precisa de los datos llegando a una precisión por encima del 99% tanto en el conjunto de prueba, como de validación

Como mencionamos anteriormente los supuestos no se cumplen por lo que no podemos asegurar que los resultados sean totalmente confiables

Por otro lado, los porcentajes de comparacion entre prueba y entrenamiento son sospechosamente altos, por lo que podriamos pensar que se esta generando un sobreajuste a pesar de separar el conjunto de datos en entretamiento y prueba

Sería adecuado poder hacer este análisis con un conjunto de datos mucho mayor, para poder tener una mayor confiabilidad en el modelo

Casualmente dan iguales o mejores resultados en el conjunto de prueba que en el de entrenamiento, considero que esto se da por la alta precisión del modelo y la baja cantidad de datos para probar

4.- Elija un método de Clasificación jerárquico y aplíquelo a los datos. Interprete los resultados.

Elijo un subset de 1000 datos para realizar el dendograma y utilizo todos los atributos del mismo para su el método de clasificación jerárquico, virifico que los subsets esten balanceados
df_clustering = sample_data
cantidad_clusters=2
set.seed(661)
par(mfcol = c(1,1))
data_c_diag = df_clustering[-c(13)] # 
sample_cluster1 <- data_c_diag[sample(1:nrow(data_c_diag), 1000,replace=FALSE),]
sample_cluster <- as.data.frame(scale(sample_cluster1[,1:12]))
sample_cluster$variedad <- sample_cluster1$variedad

# Imprimo cuántos valores hay en cada categoría de la variable diagnóstico en mi muestra de n=100
kable(table(sample_cluster$variedad))
Var1 Freq
1 503
2 497
Escogemos la distancia euclidea para realizar la matriz de distancias, y luego comparamos los diferentes métodos de distancias entre grupos para quedarnos con el más adecuado
# Escalo los datos y hago PCA
df_sample_cluster = sample_cluster[-13]
df_sample_cluster_pca <-df_sample_cluster %>% dplyr::select(is.numeric)
df_clustering.pc2 = prcomp(df_sample_cluster_pca,scale = TRUE)
# Matriz de distancias euclídeas 
mat_dist <- dist(x = df_sample_cluster, method = "euclidean") 
# Dendrogramas (según el tipo de segmentación jerárquica aplicada)  
hc_complete <- hclust(d = mat_dist, method = "complete") 
hc_average  <- hclust(d = mat_dist, method = "average")
hc_single   <- hclust(d = mat_dist, method = "single")
hc_ward     <- hclust(d = mat_dist, method = "ward.D2")
Calculo el coeficiente de correlación cofenético saber con que método avanzar
#Calculo del coeficiente de correlacion cofenetico
completo <- round(cor(x = mat_dist, cophenetic(hc_complete)),3)
promedio <- round(cor(x = mat_dist, cophenetic(hc_average)),3)
simple <- round(cor(x = mat_dist, cophenetic(hc_single)),3)
ward <- round(cor(x = mat_dist, cophenetic(hc_ward)),3)
valores_coef <- cbind(completo,promedio,simple,ward)

# Imprimo valores de coeficiente cofenético
kable(valores_coef)
completo promedio simple ward
0.587 0.772 0.649 0.431
A pesar de que el mayor índice cofenético lo arrojo el método promedio, clusterizamos con los otros tres para poder comparar los dendogramas
# Armo clusters
jer_ward<-cutree(hc_ward,k=cantidad_clusters)           
jer_average<-cutree(hc_average,k=cantidad_clusters)      
jer_complete<-cutree(hc_complete,k=cantidad_clusters)           
jer_single<-cutree(hc_single,k=cantidad_clusters)     
# Agrego cluster a dataframe
sample_cluster$jer_ward=jer_ward
sample_cluster$jer_average=jer_average
sample_cluster$jer_complete=jer_complete
sample_cluster$jer_single=jer_single
mar = c(5.1, 4.1, 4.1, 2.1) 
pch=c('royalblue2','#ff7474ff') 
cols=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_ward))]],0.7)
dend_ward <- color_branches(as.dendrogram(hc_ward), k = 2)
dend_ward <- set(dend_ward, "labels_cex", 0.1)
grafico1 <- dend_ward %>%  set("leaves_pch",19)%>%
        set("leaves_cex", .9) %>% set("leaves_col", cols) %>% 
        plot(main = "Dendrograma jerárquico", ylab='Distancia',cex.lab=1, cex.axis=.6)+
        mtext(side = 3, line = 0.5, at = 1, adj = -1.8, 'Distancia Ward')+
        mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(5,28, title='Variedad', 
     legend = c("blanco" , "tinto"), 
     col = c('royalblue2','#ff7474ff') , 
     pch = c(19,19), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0, 0.1))

cols_a=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_average))]],0.7)
dend_average <- color_branches(as.dendrogram(hc_average), k = 2)
dend_average <- set(dend_average, "labels_cex", 0.1)
grafico2 <- dend_average %>% set("leaves_pch",19) %>%  
        set("leaves_cex", .8) %>%  set("leaves_col", cols_a) %>% 
        plot(main = "Dendrograma jerárquico",  ylab='Distancia',cex.lab=1, cex.axis=.6)+
        mtext(side = 3, line = 0.5, at = 1, adj = -1.3, 'Distancia Promedio')+
        mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(76,8, title='Variedad', 
     legend = c("blanco" , "tinto"), 
     col = c('royalblue2','#ff7474ff') , 
     pch = c(19,19), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0, 0.1))

cols_c=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_complete))]],0.7)
dend_complete <- color_branches(as.dendrogram(hc_complete), k = 2)
dend_complete <- set(dend_complete, "labels_cex", 0.1)
grafico3 <- dend_complete %>% set("leaves_pch",19) %>%  # node point type
        set("leaves_cex", .8) %>%  # node point size
        set("leaves_col", cols_c) %>% #node point color
        plot(main = "Dendrograma jerárquico", ylab='Distancia',cex.lab=1, cex.axis=.6)+
        mtext(side = 3, line = 0.5, at = 1, adj = -1.3, 'Distancia Completa')+
        mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(85,16, title='Diagnóstico', 
     legend = c("blanco" , "tinto"), 
     col = c('royalblue2','#ff7474ff') , 
     pch = c(19,19), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0, 0.1))

cols_s=alpha(pch[sample_cluster$variedad[order.dendrogram(as.dendrogram(hc_single))]],0.7)
dend_single <- color_branches(as.dendrogram(hc_single), k = 2)
dend_single <- set(dend_single, "labels_cex", 0.1)
grafico4 <- dend_single %>% set("leaves_pch",19) %>%  # node point type
        set("leaves_cex", .8) %>%  # node point size
        set("leaves_col", cols_s) %>% #node point color
        plot(main = "Dendrograma jerárquico", ylab='Distancia',cex.lab=1, cex.axis=.7)+
        mtext(side = 3, line = 0.5, at = 1, adj = -1.7, 'Distancia Simple')+
        mtext(side = 1, line = 0.5, at = 1, adj = -4.1, 'Vinos')
legend(80,7, title='Variedad', 
     legend = c("blanco" , "tinto"), 
     col = c('royalblue2','#ff7474ff') , 
     pch = c(19,19), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0, 0.1))

# ·····················································
# cuántos pacientes de cada diagnóstico están en cada cluster:
ward_cluster1 <- sample_cluster %>% filter (jer_ward == '1')
cluster1 <- table(ward_cluster1$variedad)
Ward_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
ward_cluster2 <- sample_cluster %>% filter (jer_ward == '2')
cluster2 <- table(ward_cluster2$variedad)
Ward_cluster.2 <- round(prop.table(cluster2)*100,2)

kable(cbind(rbind(cluster1,cluster2),rbind(Ward_cluster.1,Ward_cluster.2)))
1 2 1 2
cluster1 32 490 6.13 93.87
cluster2 471 7 98.54 1.46 
# ·····················································
promedio_cluster1 <- sample_cluster %>% filter (jer_average == '1')
cluster1 <- table(promedio_cluster1$variedad)
promedio_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
promedio_cluster2 <- sample_cluster %>% filter (jer_average == '2')
cluster2 <- table(promedio_cluster2$variedad)
promedio_cluster.2 <- round(prop.table(cluster2)*100,2)

kable(cbind(rbind(cluster1,cluster2),rbind(promedio_cluster.1,promedio_cluster.2)))
1 2 1 2
cluster1 502 497 50.25 49.75
cluster2 1 0 100.00 0.00 
# ·····················································
completo_cluster1 <- sample_cluster %>% filter (jer_complete == '1')
cluster1 <- table(completo_cluster1$variedad)
completo_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
completo_cluster2 <- sample_cluster %>% filter (jer_complete == '2')
cluster2 <- table(completo_cluster2$variedad)
completo_cluster.2 <- round(prop.table(cluster2)*100,2)

kable(cbind(rbind(cluster1,cluster2),rbind(completo_cluster.1,completo_cluster.2)))
1 2 1 2
cluster1 502 497 50.25 49.75
cluster2 1 0 100.00 0.00 
# ·····················································
simple_cluster1 <- sample_cluster %>% filter (jer_single == '1')
cluster1 <- table(simple_cluster1$variedad)
simple_cluster.1 <- round(prop.table(cluster1)*100,2)
# ·····················································
simple_cluster2 <- sample_cluster %>% filter (jer_single == '2')
cluster2 <- table(simple_cluster2$variedad)
simple_cluster.2 <- round(prop.table(cluster2)*100,2)


kable(cbind(rbind(cluster1,cluster2),rbind(simple_cluster.1,simple_cluster.2)))
1 2 1 2
cluster1 502 497 50.25 49.75
cluster2 1 0 100.00 0.00

Conclusiones Como resultado del análisis de clasificación realizada utilizando diferentes tipos de segmentación jerárquica, determinamos que a pesar de que el cálculo del coeficiente de correlación cofenético nos arrojo con mayor valor al tipo promedio, ejecutando el denograma de cada uno vemos que claramente el que mejor realiza la separación de clusteres es el tipo ward, arrojando resultados muchísimos mejores que el resto

1 2 1 2 cluster1 32 490 6.13 93.87 cluster2 471 7 98.54 1.46

Limito el numero de clusters por el costo computacional

5.- Aplique a los datos el método de clasificación no jerárquico K-means. Interprete los resultados.

Realizamos el análisis correspondiente para elegir la cantidad de clusters a calcular con kmeans. Graficamos las métricas SIL y SSE
#datos_para_kmeans = sample_data[1:12]
df_kmeans = sample_data
datos_para_kmeans <- sample_data %>% dplyr::select(is.numeric)

#analisis de la cantidad de clusters. Este primer bloque es solo para definir funciones.
#se define una funcion para calcular metricas que orientan sobre el numero de clusters a elegir para el problema.
metrica = function(datA_esc,kmax,f) {
        sil = array()
        sse = array()
        datA_dist= dist(datA_esc,method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
        for ( i in  2:kmax) { if (strcmp(f,"kmeans")==TRUE) {   #centroide: tipico kmeans
                        CL  = kmeans(datA_esc,centers=i,nstart=50,iter.max = kmax)
                        sse[i]  = CL$tot.withinss 
                        CL_sil = silhouette(CL$cluster, datA_dist)
                        sil[i]  = summary(CL_sil)$avg.width}
                if (strcmp(f,"pam")==TRUE){       #medoide: ojo porque este metodo tarda muchisimo 
                        CL = pam(x=datA_esc, k=i, diss = F, metric = "euclidean")
                        sse[i]  = CL$objective[1] 
                        sil[i]  = CL$silinfo$avg.width}}
        sse
        sil
        return(data.frame(sse,sil))}
#en este bloque se estudia cuantos clusters convendría generar segun indicadores tipicos -> por ejemplo el "Silhouette"
kmax = 10

m1   = metrica(scale(datos_para_kmeans),kmax,"kmeans")  #tipica con estimadores de la normal
m1 <- m1[complete.cases(m1),]
m1$kcluster <- seq(2,kmax,1)
m1 <- m1%>%dplyr::select(3,1,2)
m1_sse <- m1%>%dplyr::select(-3)%>%mutate(metric='SSE')
colnames(m1_sse) <- c('kcluster','value','metric')
m1_sil <- m1%>%dplyr::select(-2)%>%mutate(metric='SIL')
colnames(m1_sil) <- c('kcluster','value','metric')
m1 <- rbind(m1_sse,m1_sil)
# Grafico de métricas SIL y SSE
ggplot(m1, aes(kcluster, value, linetype=metric)) + geom_line(col='red') + 
        facet_wrap(~metric, ncol=1, scales='free')+theme+geom_point(col='red', size=2, fill='pink', shape=21)+
        labs(title='Determinación de número de clusters', 
             x='k Número de clusters', y='Valor', linetype='Métrica')+
        scale_x_continuous(breaks = seq(1, kmax, by = 1))+
        scale_linetype_manual(values=c(1,2))

A sabiendas que nuestra variable objetivo tiene solo dos categorías, realizamos la clusterización con kmeans = 2, y graficamos los resultados en relación a las variables pH, densidad y calidad.
set.seed(661)
cantidad_clusters=2

CL  = kmeans(scale(datos_para_kmeans),cantidad_clusters)
df_kmeans$kmeans = CL$cluster

# Grafico scatterplot original + cluster con k=2
par(mfrow=c(1,2))
# -------------------------------------------------------------------------------------------
col1 <- c('royalblue2','#ff7474ff')
col1 <- col1[as.numeric(df_kmeans$variedad)]

scatterplot3d(df_kmeans$pH,df_kmeans$densidad,df_kmeans$calidad, color = alpha(col1,0.3), box=F,angle=45, pch = 19, grid = TRUE, tick.marks = FALSE, xlab = "ph", ylab = "densiidad", zlab = "calidad", main='Realidad')
#legend("topright", bty = "n", cex = .9, title = "Variedad", c("blanco", "tinto"), fill = c('royalblue2','#ff7474ff'))
# -------------------------------------------------------------------------------------------
colors <- c('orange','#a25da2a5')
colors <- colors[as.numeric(df_kmeans$kmeans)]

scatterplot3d(df_kmeans$pH,df_kmeans$densidad,df_kmeans$calidad, color = alpha(colors,0.3), box=F,angle=45, pch = 19, grid = TRUE, tick.marks = FALSE, xlab = "ph", ylab = "densiidad", zlab = "calidad", main='Clustering')

#legend("topright", bty = "n", cex = .9, title = "Grupo k-means", c("1", "2"), fill = c('orange','#a25da2a5'))
#conviene en un biplot ya que tengo las flechas de las variables originales
# GRAFICO ORIGINAL
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1,alpha=0.5,
         groups=factor(df_kmeans$variedad)) +
        scale_color_manual(name="variedad",values=c('royalblue2','#ff7474ff'),
                           labels=c("blanco",'tinto')) +
        theme + labs(title='Análisis de componentes principales') + 
        theme(legend.position=c(.85,.15)) + 
        labs(x='PC1 (51.5% explicado)', y= 'PC2 (22.5% explicado)')

# -------------------------------------------------------------------------------------------
# GRAFICO KMEANS
ggbiplot(df_svm.pc, obs.scale=0.1 ,var.scale=1, alpha=0.5,groups = as.factor(df_kmeans$kmeans) )+
        scale_color_manual(name="Grupo k-means", values=c("orange",'#a25da2a5',"darkgreen"),
                           labels=c("1", "2","3")) + theme+
        labs(title= 'Representación del clustering utilizando k-means')+theme(legend.position=c(.85,.15)) 

# ·····················································
# cuántos pacientes de cada diagnóstico están en cada cluster:
pacientes_cluster1 <- df_kmeans %>% filter (kmeans == '1')
cluster1 <- table(pacientes_cluster1$variedad)
porcentaje_cluster.1 <- round(prop.table(cluster1)*100,2)
pacientes_cluster2 <- df_kmeans %>% filter (kmeans == '2')
cluster2 <- table(pacientes_cluster2$variedad)
porcentaje_cluster.2 <- round(prop.table(cluster2)*100,2)
# ·····················································
# Imprimo resultados
kable(cbind(rbind(cluster1,cluster2),rbind(porcentaje_cluster.1,porcentaje_cluster.2)))
1 2 1 2
cluster1 980 13 98.69 1.31
cluster2 20 987 1.99 98.01 
kmeans1 <- df_kmeans %>% filter(kmeans==1)  %>%dplyr::select(c(1:12))%>% colMeans()
kmeans2 <- df_kmeans %>% filter(kmeans==2)  %>%dplyr::select(c(1:12))%>% colMeans()
blanco <- df_kmeans %>%filter(variedad=='1') %>%dplyr::select(c(1:12))%>% colMeans()
tinto <- df_kmeans %>%filter(variedad=='2') %>%dplyr::select(c(1:12))%>% colMeans()
# Imprimo resultados
kable(rbind(kmeans1,blanco,kmeans2,tinto))
acidez_fija acidez_volatil acido_citrico azucar_residual cloruros anhidrido_sulfuroso_libre anhidrido_sulfuroso_total densidad pH sulfatos alcohol calidad
kmeans1 6.857150 0.2776788 0.3430111 6.538872 0.0459849 36.41138 141.04632 0.9940956 3.182699 0.4889325 10.49300 5.905337
blanco 6.853950 0.2820350 0.3403700 6.490900 0.0462090 35.94800 140.12700 0.9940928 3.185840 0.4907700 10.48465 5.871000
kmeans2 8.340119 0.5239623 0.2740516 2.529146 0.0876207 15.59533 46.08937 0.9967512 3.311668 0.6625025 10.44361 5.612711
tinto 8.353700 0.5213300 0.2762100 2.549050 0.0876880 15.91300 46.34400 0.9967726 3.309430 0.6618800 10.45162 5.645000 
# ya sabíamos del qqplot que las variables no son normales univariadas (si se separan por diagnóstico). por TLC, no obstante, uno puede aproximar las univariadas y hacer un test acorde para ver diferencias de medias

blanco_N <- df_kmeans%>%filter(variedad=='1')
tinto_M <- df_kmeans%>%filter(variedad=='2')

pH <- z.test(blanco_N$pH, sigma.x=sd(blanco_N$pH), tinto_M$pH, sigma.y=sd(tinto_M$pH), conf.level=0.95)$p.value
densidad <- z.test(blanco_N$densidad, sigma.x=sd(blanco_N$densidad), tinto_M$densidad, sigma.y=sd(tinto_M$densidad), conf.level=0.95)$p.value
calidad <- z.test(blanco_N$calidad, sigma.x=sd(blanco_N$calidad), tinto_M$calidad, sigma.y=sd(tinto_M$calidad), conf.level=0.95)$p.value

blanco_1 <- df_kmeans%>%filter(kmeans=='1')
tinto_2 <- df_kmeans%>%filter(kmeans=='2')

pH_c <- z.test(blanco_1$pH, sigma.x=sd(blanco_1$pH), tinto_2$pH, sigma.y=sd(tinto_2$pH), conf.level=0.95)$p.value
densidad_c <- z.test(blanco_1$densidad, sigma.x=sd(blanco_1$densidad), tinto_2$densidad, sigma.y=sd(tinto_2$densidad), conf.level=0.95)$p.value
calidad_c <- z.test(blanco_1$calidad, sigma.x=sd(blanco_1$calidad), tinto_2$calidad, sigma.y=sd(tinto_2$calidad), conf.level=0.95)$p.value

variable <- c('variedad','cluster')
# Imprimo resultados
kable(rbind(variable,cbind(rbind(pH, densidad, calidad),rbind(pH_c, densidad_c, calidad_c))))
variable variedad cluster
pH 7.24729192241589e-73 2.36190019262579e-80
densidad 5.95484451106308e-126 1.94113514456475e-122
calidad 5.17424960452661e-09 2.88274220391404e-14 
#datos_para_kmeans = sample_data[1:12]
df_kmeans_4 = sample_data
datos_para_kmeans_4 <- sample_data %>% dplyr::select(is.numeric)

# ENTRENO K MEANS CON K=3
cantidad_clusters_4=4
set.seed(661)
CL_4  = kmeans(scale(datos_para_kmeans_4),cantidad_clusters_4)
df_kmeans_4$kmeans = CL_4$cluster
colors2 <- c('orange','#a25da2a5', 'darkgreen','lightblue')
colors2 <- colors2[as.numeric(df_kmeans_4$kmeans)]
# -------------------------------------------------------------------------------------------
# GRAFICO K=3
scatterplot3d(df_kmeans_4$pH,df_kmeans_4$densidad,df_kmeans_4$calidad, color = alpha(colors2,0.3), box=F,angle=25, pch = 19, grid = TRUE, tick.marks = FALSE, xlab = "pH", ylab = "densidad", zlab = "calidad", main='Clustering k=4', cex.lab=.8,scale.y=2, cex.symbols=1.3)
legend(x=4, y=5.5, bty = "n", cex = 1.1, title = "Grupo k-means", c("1", "2","3","4"), fill = c('orange','#a25da2a5', 'darkgreen','lightblue'))

# ·····················································
# cuántos pacientes de cada diagnóstico están en cada cluster:
pacientes_cluster1_4 <- df_kmeans_4 %>% filter (kmeans == '1')
cluster1_4 <- table(pacientes_cluster1_4$variedad)
porcentaje_cluster_4.1 <- round(prop.table(cluster1_4)*100,2)
pacientes_cluster2_4 <- df_kmeans_4 %>% filter (kmeans == '2')
cluster2_4 <- table(pacientes_cluster2_4$variedad)
porcentaje_cluster_4.2 <- round(prop.table(cluster2_4)*100,2)

pacientes_cluster3_4 <- df_kmeans_4 %>% filter (kmeans == '3')
cluster3_4 <- table(pacientes_cluster3_4$variedad)
porcentaje_cluster_4.3 <- round(prop.table(cluster3_4)*100,2)
pacientes_cluster4_4 <- df_kmeans_4 %>% filter (kmeans == '4')
cluster4_4 <- table(pacientes_cluster4_4$variedad)
porcentaje_cluster_4.4 <- round(prop.table(cluster4_4)*100,2)

# ·····················································
# Imprimo resultados
kable(cbind(rbind(cluster1_4,cluster2_4,cluster3_4,cluster4_4),rbind(porcentaje_cluster_4.1,porcentaje_cluster_4.2,porcentaje_cluster_4.3,porcentaje_cluster_4.4)))
1 2 1 2
cluster1_4 415 3 99.28 0.72
cluster2_4 15 548 2.66 97.34
cluster3_4 9 397 2.22 97.78
cluster4_4 561 52 91.52 8.48 
kmeans_41 <- df_kmeans_4 %>% filter(kmeans==1)  %>%dplyr::select(c(1:12))%>% colMeans()
kmeans_42 <- df_kmeans_4 %>% filter(kmeans==2)  %>%dplyr::select(c(1:12))%>% colMeans()
kmeans_43 <- df_kmeans_4 %>% filter(kmeans==3)  %>%dplyr::select(c(1:12))%>% colMeans()
kmeans_44 <- df_kmeans_4 %>% filter(kmeans==4)  %>%dplyr::select(c(1:12))%>% colMeans()

blanco <- df_kmeans_4 %>%filter(variedad=='1') %>%dplyr::select(c(1:12))%>% colMeans()
tinto <- df_kmeans_4 %>%filter(variedad=='2') %>%dplyr::select(c(1:12))%>% colMeans()
# Imprimo resultados
kable(rbind(kmeans_41,kmeans_42,kmeans_43,kmeans_44,blanco,tinto))
acidez_fija acidez_volatil acido_citrico azucar_residual cloruros anhidrido_sulfuroso_libre anhidrido_sulfuroso_total densidad pH sulfatos alcohol calidad
kmeans_41 6.927990 0.2887560 0.3575598 11.153469 0.0539545 45.71651 169.04426 0.9969552 3.157679 0.4869378 9.531738 5.543062
kmeans_42 7.378686 0.6141297 0.1364121 2.399290 0.0808490 16.28686 49.15631 0.9963379 3.376927 0.5980462 10.196773 5.358792
kmeans_43 9.816749 0.4133498 0.4635468 2.729926 0.0994852 14.54680 42.41626 0.9976136 3.208054 0.7437192 10.567406 5.884237
kmeans_44 6.805791 0.2758401 0.3297227 3.129935 0.0414927 28.83524 115.68434 0.9921187 3.216444 0.5064600 11.290131 6.187602
blanco 6.853950 0.2820350 0.3403700 6.490900 0.0462090 35.94800 140.12700 0.9940928 3.185840 0.4907700 10.484650 5.871000
tinto 8.353700 0.5213300 0.2762100 2.549050 0.0876880 15.91300 46.34400 0.9967726 3.309430 0.6618800 10.451617 5.645000

Conclusiones Ejecutamos el modelo k means con el conjunto de datos, para esto en primer lugar realizamos el análisis correspondiente para elegir la cantidad de clusters a calcular. para esto calculamos las métricas Silhouette y SSE, obtenemos que el conjunto de datos seria bueno evaluarlo con 4, 5 clusters, sabiendo que nuestra variable target a dividir tan solo tiene dos categorías, continuamos la implementación con k means = 2 y k means = 4 para comparar resultados. Graficando un biplot con lo resuelto por K = 2 vemos que genera dos grupos muy marcados, donde no hay prácticamente un solapamiento entre ambos lo que sí sucede conjunto original, por el contrario con K = 4, si vemos que las fronteras entre los grupos comienzan a solaparse, por lo que con un análisis más profundo, la agrupación de varios clusters a través del análisis de sus medias podría brindar una imagen más acorde a la realidad. En este caso, observamos que k means = 1 y k means = 4 tienen medias muy similares, que podríamos asociarlas a la distribución de vino blanco.

C.- Presente un informe final de 2 carillas como máximo, no incluya gráficos, explicando las conclusiones del trabajo realizado, mencione si es necesario validar supuestos requeridos para aplicar el método. Compare los resultados de los métodos supervisados y establezca conclusiones. Por otro lado, compare los métodos no supervisados y presente sus conclusiones.

D.- Realice el trabajo en R y adjunte el código en R, el RMD y el html donde se puede observar el código y las salidas junto con la interpretación de cada método.